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Typical quasistatic compression algorithms for generating jammed packings of athermal, purely 
repulsive particles begin with dilute configurations and then apply successive compressions with 
relaxation of the elastic energy allowed between each compression step. It is well-known that during 
isotropic compression athermal systems with purely repulsive interactions undergo a jamming tran- 
sition at packing fraction <f>j from an unjammed state with zero pressure to a jammed, rigid state 
with nonzero pressure. Using extensive computer simulations, we show that a novel second-order- 
like transition, the contact percolation transition, which signals the formation of a system-spanning 
cluster of mutually contacting particles, occurs at <f>p < </>j, preceding the jamming transition. By 
measuring the number of non-floppy modes of the dynamical matrix, and the displacement field and 
time-dependent pressure following compression, we find that the contact percolation transition also 
heralds the onset of complex spatiotemporal response to applied stress. Thus, highly heterogeneous, 
cooperative, and non-affine particle motion occurs in unjammed systems significantly below the 
jamming transition for 4>p < (j> < <f>j, not only for jammed systems with (j> > 4>.j. 
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The jamming transition in athermal, purely repulsive 
particulate systems, such as granular media, foams, and 
colloids, has been characterized extensively in computer 
simulations [l[ and experiments 0. For example, when 
model frictionless spheres are compressed to packing frac- 
tions above <f>j, static particle configurations transition 
from unjammed to jammed with nonzero pressure (and 
elastic energy) and a rigid backbone of contacting par- 
ticles with all positive nontrivial eigenvalues Q of the 
dynamical matrix. Signatures of the jamming transi- 
tion, such as anomalous scaling of the zero-frequency 
shear modulus with packing fraction Q and diverging 
length scales || associated with cooperative particle 
rearrangements, have been investigated in thermal sys- 
tems |7| in the zero-temperature limit and in sheared 
systems [8l4llj in the zero-shear rate limit, but mainly 
for packing fractions near and above <f>j. 

However, there have been few detailed studies of the 
structural and mechanical properties of unjammed ather- 
mal particulate systems well below (j)j. As shown in 
Fig. [TJ typical quasistatic compression algorithms used to 
generate static packings in experiments start with a di- 
lute collection of particles, and the sample is successively 
compressed by small amounts with energy relaxation al- 
lowed between each compression step. For <j> < tfrj, the 
configurations are not completely rigid, and thus after 
each small compression, particles can rearrange until all 
interparticle forces are zero. Despite this, the response 
of unjammed packings to compression and other pertur- 
bations can be highly heterogeneous, cooperative, and 
non-affine. Thus, an important, unanswered question is 
at what packing fraction does complex spatiotemporal 
response first occur in unjammed particulate systems? 

In this letter, we describe computational studies of a 



novel second-order-like transition — the contact percola- 
tion transition at 4>p — in athermal particulate systems of 
purely repulsive, frictionless disks that signals the forma- 
tion of a system-spanning cluster of connected non-force- 
bearing interparticle contacts and the onset of complex 
spatiotemporal response well below the jamming transi- 
tion. These systems display robust scaling behavior near 
4>p, but with a correlation length exponent that differs 
from the corresponding values for random continuum [l2j 
and rigidity percolation fl3j ] . In addition, we find that 
the the number of 'blocked' degrees of freedom, and the 
particle displacement and stress relaxation time in re- 
sponse to compression begin to increase significantly for 
4> > 4>p. These results emphasize that highly coopera- 
tive and glassy dynamics occur in compressed athermal 
systems significantly below the jamming transition, not 
only for jammed systems with <f> > <f>j . 

We focus on systems composed of bidisperse friction- 
less disks that interact via the purely repulsive linear 
spring potential: 



(1) 



where e is the characteristic energy scale, 



is the 



Heaviside step function, m is the separation between the 
centers of disks i and j , and cry = (o^ + aj ) /2 is their av- 
erage diameter. We chose a 50 — 50 (by number) mixture 
of large and small disks with diameter ratio ai/a s = 1.4 
to inhibit crystallization [lj] and implemented periodic 
boundary conditions in a unit square. We employed a 
quasistatic isotropic compression algorithm to gen erate 
static packings over a range of packing fractions [l4[ . We 
initialize each system with random particle positions at 
<\> = and zero velocities. We then compress the sys- 
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FIG. 1: (Color online) Typical snapshots from the quasistatic isotropic compression algorithm to generate static particle 
configurations as a function of packing fraction <f>. Particles with a given color shading belong to the same cluster of mutually 
contacting particles. For (left) <j> < <f>p, the system is unjammed, and the largest cluster of contacting particles does not 
percolate. The largest cluster begins to percolate for (middle) cj>p < cj> < <f)j, but the system remains unjammed since it 
possesses nontrivial floppy modes. For (right) <j> > <j>j, the percolating cluster of contacting particles has become rigid except 
for a small number of rattlers and the system is jammed. 



tem in steps of Acj> = 1CP 3 and relax the small particle 
overlaps after each step by solving Newton's equations of 
motion in the overdamped limit, 



ma; 
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(2) 



where m and Hi are the particle mass and acceleration, 
F(rij) = —dV(rij)/drijfij, fij is the unit vector connect- 
ing the centers of particles i and j, and b = ba s / \/me is 
the damping coefficient, until the total potential energy 
per particle falls below a specified tolerance V/eN < 
Vtoi = 10~ 16 . We studied two values for the damping 
coefficient, 6=1 and b — > oo using steepest descent 
dynamics. We continue compression steps followed by 
relaxation until the systems jam at a configuration de- 
pendent 4>j rj 0.84. 

In Fig.f2J we characterize the contact percolation tran- 
sition by plotting the probability P(4>) that the system 
forms a system-spanning network of interparticle con- 
tacts in either the x- or y-direction, where contact is de- 
termined by Tij < cry, at each </> immediately following 
a compression step. We find that the contact percola- 
tion transition at <j)p — 0.549 <C 4>j becomes sharper 
with system size and obeys finite-size scaling, but with 
a correlation length exponent v ss 1.68 that is signifi- 
cantly larger than that for random continuum 12j and 
rigidity percolation jl3|, but smaller than that found for 
contact percolation in athermal particulate systems with 
short-range attractions [HI . (Note that percolation onset 
occurs at a similar value 4>p = 0.558 ± 0.008 for ather- 
mal systems with short-range attractions.) In contrast, 
the exponent r « 2.01 that characterizes the power-law 
scaling of the cluster size distribution and the fractal di- 



mension D sa 1.89 are similar to that for random con- 
tinuum percolation and contact percolation for athermal 
systems with short-range attractions, and obey hyper- 
scaling D(t - 1) = 2. (See Table H) We find that these 
results are insensitive to the damping coefficient b in the 
overdamped limit and compression step for A</> < 10~ 3 . 

We have shown that immediately following each com- 
pression step, a system-spanning cluster of interparticle 
contacts forms at t = for <j) > <pp 1 much below <pj. 
To determine if this purely geometrical transition influ- 
ences the mechanical properties of the system, we mea- 
sure the eigenvalues of the dynamical matrix at t = 
following each compression step. As static packings are 
compressed, they progressively become less floppy, i.e. 
fewer single and collective particle motions cost zero en- 
ergy. We quantify the increase in rigidity by measuring 
the fraction of non-floppy or 'blocked' eigenmodes F(</)) — 
the ratio of the number N n f of non-zero eigenvalues of 
the dynamical matrix to the total number of nontrivial 
modes 2N' — 2 (where N' — N — N r and N r is the num- 
ber of rattler particles at jamming [lfjl ]) — over a range of 
packing fractions. With this definition, F(<fij) — 1. In 
Fig. [3] (a) , we show that the fraction of non- floppy modes 
F((f>) grows linearly with cj) for small <p and near jamming 
F(4>) scales as 1 - A(<j)j - 0) Q with A s» 0.99 and a m 0.6. 
We find that F(<j>) begins to deviate from linear behav- 
ior near <pp. In addition, near <pp, F{(f>) does not obey 
the scaling behavior near <pj. Thus, we have identified a 
signature in the mechanical response near (pp. 

In addition, we find that the contact percolation tran- 
sition at <pp signals the onset of complex spatiotempo- 
ral response to isotropic compression. We measure two 
quantities that characterize cooperative particle motion: 
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Percolation type 



D 



Repulsive contact 


1.70 


± 


0.09 


2.01 


± 


0.04 


1.8 


9± 


0.03 


Attractive contact 


1.92 


± 


0.03 


2.04 


± 


0.04 


1.8 


8± 


0.04 


Continuum 


1.34 


± 


0.02 


2.02 


± 


0.03 


1.91 ± 


0.04 



TABLE I: Critical exponents for contact percolation in ather- 
mal systems with purely repulsive interactions and short- 
range attractive interactions [TBI ], as well as random contin- 
uum percolation [12] in 2D. 

1) the displacement in configuration space between sys- 
tems at successive compressions, 
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Xi (oo)) 2 + ( yi (Q) - yi (oo))% (3) 



where (xi(0),yi(0)) and (2^(00), yi(oo)) are the locations 
of particle i at t = and t — > 00, respectively, during 
the energy relaxation, and 2) the accumulated distance 
traveled in configuration space, 



L = I dt 
'0 



\ 



N 



(4) 



from t = after the compression to the end of the en- 
ergy relaxation. Both D and L are normalized by the 
compression step size A0. In Fig. [3] (b), we show that D 
and L grow roughly exponentially for small 0, but begin 
to deviate from the low-0 behavior near <j)p. An even 
larger signature of collective, non-affine motion is shown 
in the inset to Fig. [3] (b). The difference L — D begins 
to increase dramatically at (j>i « 0.78, which is also sig- 
nificantly below <j)j. Thus, for <f> > (f>p particles move 
collectively in response to compression, but for <j> > (j>i 
after relaxation they end up close to where they started. 

We also characterize the time-dependent response to 
compression. In Fig. |4](a), we plot the normalized pres- 
sure p = p(t)/p(0) as a function of time following a com- 
pression step A</> = 1CP 3 over a wide range of packing 
fractions. The normalized pressure develops a plateau 
near p = 0.6 and decays more slowly with increasing 
packing fraction. The plateau in the normalized pres- 
sure suggests that significant decay of the pressure in 
the system requires first the collective motion of exterior 
particles in percolating clusters followed by motion of the 
interior particles. In Fig. @] (b), we plot the pressure re- 
laxation time t* , i.e. the time at which the normalized 
pressure has decayed to a specified small value p = 10~ 4 . 
The shape of the increase in the pressure relaxation time 
is relatively insensitive to b (for overdamped dynamics) 
and p. In the inset to Fig. 2] (b), the logarithmic slope 
of t* with respect to (f> shows a subtle increase near <pp, 
but a much more dramatic increase near 4>i, where the 
the measure of cooperative motion increased significantly. 
(See the inset to Fig. [3] (b).) Thus growth in the stress 
relaxation time occurs at packing fractions significantly 
below <j)j. 




FIG. 2: Percolation probability P(<f>) that the system pos- 
sesses a system-spanning cluster (in either the x- or y- 
direction) immediately following a compression step A<f> = 
10" 3 versus packing fraction <j> for N = 100, 200, 400, 800, 
1600, 3200, and 6400 particles (from left to right) averaged 
over 400 configurations. Inset: Same as the main figure ex- 
cept the horizontal axis is scaled by ((f) — (f>p)N 1 ^ 2 " , where 
<t> P = 0.549 ± 0.001 and v = 1.70 ± 0.09. 
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FIG. 3: (a) The fraction F((f>) of non-floppy eigenmodes of the 
dynamical matrix in the system measured immediately after 
a compression step of A(f> = 10 -3 over a range of packing 
fractions for N = 1000. Fits to the low and high <j> behav- 
ior of F(<f>) are shown as dashed and solid lines, respectively, 
(b) The accumulated distance L and displacement D between 
successive compressions normalized by A(f> as a function of 
packing fraction. The \ow-(j>, roughly exponential behavior 
of L and D is indicated by solid and dashed lines. The in- 
set shows L — D versus (f>. In both panels, the vertical line 
indicates the percolation transition (j>P = 0.549. 
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FIG. 4: (a) Normalized pressure p = p(t)/p(0) versus t fol- 
lowing a compression step Acf> = 10~ 3 with 6=1 over a range 
of packing fractions including (j> = 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 
0.65, 0.7, 0.73, 0.75, 0.78, 0.8, 0.81, and 0.82 (from left to 
right) for N = 200. The horizontal line p = 10 -4 indicates 
the normalized pressure at which the pressure relaxation time 
t* was measured, (b) The pressure relaxation timescale t*(</>) 
for systems with b — > oo (triangles) and 1 (squares) . The inset 
shows the logarithmic derivative of t* with respect to <f). 

In this letter, we described extensive computational 
studies of the geometrical and mechanical properties 
of unjammed, athermal particulate systems undergoing 
quasistatic isotropic compression. A novel aspect of 



our work is that we focused on unjammed rather than 
jammed configurations (T^]. We first showed that these 
systems undergo the 'contact percolation' transition at 
fap < <t>j, much below the jamming transition. Above 
4>p, these systems possess a system-spanning cluster of 
non-force bearing interparticle contacts. Near the transi- 
tion, the cluster size distribution displays robust scaling 
with system size with a correlation length exponent v 
that is larger [l5| than that for random continuum and 
rigidity percolation, and thus contact percolation belongs 
to a distinct universality class. We also performed mea- 
surements to assess the mechanical properties and time- 
dependent stress response for unjammed systems with 
fap < 4> < faj. We found three key results: 1) the in- 
crease in the fraction F of non-floppy modes begins to 
deviate from the low-</) behavior near <pp and the ac- 
celeration of non-floppy modes near 4>p is distinct from 
the approach to F — 1 near <pj. 2) The </>-dependcncc 
of two measures of cooperative motion — the accumulated 
distance L and displacement D in configuration space be- 
tween successive compressions — deviates from the low-<^> 
behavior near (pp. Further, their difference L — D begins 
to grow dramatically near fa w 0.78, which signals block- 
ing in configuration space well below <f>j. 3) The growth 
in the stress relaxation time t* also begins to deviate from 
the low-</> behavior near fa. Thus, the process of rigidifi- 
cation and the onset of cooperative, heterogeneous, and 
slow dynamics occurs well below jamming onset in ather- 
mal particulate systems. 
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